library(ggplot2)
library(scales)
library(apsrtable)

setwd("/Users/nilsw/Projects/EventDataBias/R")

# read data
pairs <- read.csv("pairs_replication.csv", header=T)


# plot spatial error distribution (Fig. 2)
ggplot(pairs, aes(pairs$distance/1000)) + stat_ecdf(size=1.5) + xlim(0,200) + theme_bw() + xlab("Spatial error (km)") + ylab("Proportion")
ggsave("errordist.pdf", width=7, height=5)

# plotting the magnitude of the bias (Fig. 3)
pairs$ucdp_bias <- ifelse(pairs$totalkia<pairs$low,pairs$totalkia-pairs$low,0)
pairs$ucdp_bias <- ifelse(pairs$totalkia>pairs$high,pairs$totalkia-pairs$high,pairs$ucdp_bias)
ggplot(pairs[abs(pairs$ucdp_bias)<=25,]) + geom_histogram(aes(x=ucdp_bias, y=..density..), binwidth=3, origin=-25.5) + xlab("") + ylab("Proportion") + ylim(0,0.4) + annotate("text", label="UCDP = SIGACTS\n(N=544)", x=0, y=0.35) + annotate("text", label="UCDP > SIGACTS\n(N=348)", x=-15, y=0.25) + annotate("text", label="UCDP < SIGACTS\n(N=188)", x=15, y=0.25) + theme_bw()
ggsave("casualty_bias.pdf", width=7, height=5)

# Regressions
pairs$type <- relevel(as.factor(pairs$type), "Explosive Hazard")

# how location affects spatial accuracy (Table 2)
m1 <- glm(log(distance) ~ log(dist_nearest_town) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs)     
m2 <- glm(log(distance) ~ log(dist_nearest_city) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs)   
m3 <- glm(log(distance) ~ log(wd_pop+1) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs)   
apsrtable(m1, m2, m3, stars="default")


# how location affects casualty accuracy (Table 3)
m4 <- glm(incorrectcasualties ~ log(dist_nearest_town) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs, family=binomial)     
m5 <- glm(incorrectcasualties ~ log(dist_nearest_city) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs, family=binomial) 
m6 <- glm(incorrectcasualties ~ log(wd_pop+1) + as.factor(type=="Explosive Hazard") + as.factor(wd_provid), data=pairs, family=binomial) 
apsrtable(m4, m5, m6, stars="default")

